close all; clear;

% Space-time spectrum estimation from CHAMP EEJ/EEF data
% The spatio-temporal spectra is derived following the 
% paper by Wu, Dong L., Paul B.
% Hays, Wilbert R. Skinner, 1995: A Least Squares Method for Spectral
% Analysis of Space-Time Series. J. Atmos. Sci., 52, 35013511
% Instead of LS, I use robustfit function of MATLAB
% 
load /Users/manojnair/projects/tides/eef_data_2000_2009.mat eef;

% # Field 1: timestamp (UT)
% # Field 2: longitude (degrees)
% # Field 3: latitude (degrees)
% # Field 4: local time (hours)
% # Field 5: season (day of year)
% # Field 6: eastward electric field (mV/m)
% # Field 7: equatorial vertical electric field at 105 km altitude (mV/m)
% # Field 8: equatorial (UxB)_z 105 km altitude (mV/m)
% # Field 9: DC current shift (A/m)
% # Field 10: CHAMP peak current value (A/m)
% # Field 11: KP
% # Field 12: F10.7 (W/m^2)
% # Field 13: F10.7A (W/m^2)
% # Field 14: R^2 (coefficient of determination)
% # Field 15: R (correlation of CHAMP and model profiles)
% # Field 16: chi^2
% # Field 17: EEFM modeled eastward electric field (mV/m)

% find day of the year for seasonal average

datemat = datevec(eef(:,1)); %getting date vector form date number
doyvec = dayofyear(datemat(:,1)',datemat(:,2)',datemat(:,3)'); %get DOY


% From hermann Luehr : Here are some suggestions.
%  -  Please use the EEJ peak current density values (determined by Stefan) directly, not the EEF.
%  -  Provide separate spectra for the seasons (Mai, June, July, Aug.), (Nov. Dec. Jan. Feb.),
%     (15 June - 15 Oct). The last period is centered on DE3 maximum. 4 months are in any case preferable because of the good local time coverage of CHAMP.
%  -  Is your method able to provide phase information?
%  For the full spectra you use all CHAMP years. I would suggest to separate it into the first 5 years, 
% Aug. 2000-2005 and Aug. 2005-2010, 

days_in_prev_months = [0 31 59 90 120 151 181 212 243 273 304 334 366];
lt_min = 7;
lt_max = 17;
date_min = datenum(2000,08,1);
date_max = datenum(2005,07,31);
% date_min = datenum(2005,08,1);
% date_max = max(eef(:,1));

p_test_lim = 0.05;
periods = 3:1:30;
wavenum = -10:1:10;

% normalized EEJ strength
% for normalizing the effect of conductvity, I am dividing the 
% eej data with sqrt(cos((pi/12).*(LT - 12.5))))

eef(:,10) = eef(:,10)./sqrt(cos((pi/12).*(eef(:,4) - 12.5)));

% delete the complex numbers

L = (eef(:,10) == real(eef(:,10)));

eef(~L,10) = 0;

%% Monthly estimates of tides
for months = 1:12,

 st = days_in_prev_months(months);
 en = days_in_prev_months(months+1);



LL = doyvec >= st & doyvec <= en & eef(:,4)' >= lt_min & eef(:,4)' <= lt_max & eef(:,1)' > date_min & eef(:,1)' < date_max;
[n_data,n_days]=hist(doyvec(LL));

% make normalized latitudes between 0 (= 0) and 360 (=1)

if exist('lonnorm'),
    clear lonnorm;
end;

selected_lon = eef(LL,2);
L = selected_lon < 0;
lonnorm(L) = 360 + selected_lon(L) ;
lonnorm(~L) =  selected_lon(~L) ;
lonnorm = lonnorm/360;


clear spectra_ls spectra_rob_eej spectra_ls_eej p_st std_err stats
for i = 1:length(periods),
    for j = 1 : length(wavenum),
    
prd = periods(i);

x_m2 = [ cos(2*pi*((eef(LL,1)/( prd/24)) + (lonnorm'*wavenum(j)))) ...
         sin(2*pi*((eef(LL,1)/( prd/24)) + (lonnorm'*wavenum(j))))];  
[spectra_rob_eej(i,j,:), stats(i,j) ] = robustfit(x_m2,eef(LL,10));
p_st(i,j,:) = stats(i,j).p;
std_err(i,j,:) = stats(i,j).se;
%spectra_ls_eej(i,j,:) = x_m2\(eef(LL,10));
end;
end;
eval(['save /Users/manojnair/projects/tides/plots/' sprintf('DOY_%d_%d_LT_%d_%d_T_%s_%s',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29)) '_ROB_SPECTRA' ...
    ' spectra_rob_eej p_st std_err st en periods wavenum n_data n_days']);
end;

%% Estimate the tides for all days

st = 1;
en = 366;

LL = doyvec >= st & doyvec <= en & eef(:,4)' >= lt_min & eef(:,4)' <= lt_max & eef(:,1)' > date_min & eef(:,1)' < date_max;
[n_data,n_days]=hist(doyvec(LL));

% make normalized latitudes between 0 (= 0) and 360 (=1)

if exist('lonnorm'),
    clear lonnorm;
end;

selected_lon = eef(LL,2);
L = selected_lon < 0;
lonnorm(L) = 360 + selected_lon(L) ;
lonnorm(~L) =  selected_lon(~L) ;
lonnorm = lonnorm/360;



clear spectra_ls spectra_rob_eej spectra_ls_eej p_st std_err stats
for i = 1:length(periods),
    for j = 1 : length(wavenum),
    
prd = periods(i);

x_m2 = [ cos(2*pi*((eef(LL,1)/( prd/24)) + (lonnorm'*wavenum(j)))) ...
         sin(2*pi*((eef(LL,1)/( prd/24)) + (lonnorm'*wavenum(j))))];  
[spectra_rob_eej(i,j,:), stats(i,j) ] = robustfit(x_m2,eef(LL,10));
p_st(i,j,:) = stats(i,j).p;
std_err(i,j,:) = stats(i,j).se;
%spectra_ls_eej(i,j,:) = x_m2\(eef(LL,10));
end;
end;
eval(['save /Users/manojnair/projects/tides/plots/' sprintf('DOY_%d_%d_LT_%d_%d_T_%s_%s',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29)) '_ROB_SPECTRA' ...
    ' spectra_rob_eej p_st std_err st en periods wavenum n_data n_days']);
%% Estimate for  seasons

season_st = [dayofyear(2000,05,01) dayofyear(2000,06,15) dayofyear(2000,11,01)];
season_en = [dayofyear(2000,08,30) dayofyear(2000,10,15) dayofyear(2000,02,29)];

for seasons = 1:3,

 st = season_st(seasons);
 en = season_en(seasons);


 if en > st,
     LL = doyvec >= st & doyvec <= en & eef(:,4)' >= lt_min & eef(:,4)' <= lt_max & eef(:,1)' > date_min & eef(:,1)' < date_max;
 else, %seasons across year borders
     LL =  ( doyvec >= st & doyvec < 366 ) | ( doyvec > 0 & doyvec <= en ) & eef(:,4)' >= lt_min & eef(:,4)' <= lt_max & eef(:,1)' > date_min & eef(:,1)' < date_max; % for seasons across year border
 end;
[n_data,n_days]=hist(doyvec(LL));

% make normalized latitudes between 0 (= 0) and 360 (=1)

if exist('lonnorm'),
    clear lonnorm;
end;

selected_lon = eef(LL,2);
L = selected_lon < 0;
lonnorm(L) = 360 + selected_lon(L) ;
lonnorm(~L) =  selected_lon(~L) ;
lonnorm = lonnorm/360;



clear spectra_ls spectra_rob_eej spectra_ls_eej p_st std_err stats
for i = 1:length(periods),
    for j = 1 : length(wavenum),
    
prd = periods(i);

x_m2 = [ cos(2*pi*((eef(LL,1)/( prd/24)) + (lonnorm'*wavenum(j)))) ...
         sin(2*pi*((eef(LL,1)/( prd/24)) + (lonnorm'*wavenum(j))))];  
[spectra_rob_eej(i,j,:), stats(i,j) ] = robustfit(x_m2,eef(LL,10));
p_st(i,j,:) = stats(i,j).p;
std_err(i,j,:) = stats(i,j).se;
end;
end;
eval(['save /Users/manojnair/projects/tides/plots/' sprintf('DOY_%d_%d_LT_%d_%d_T_%s_%s',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29)) '_ROB_SPECTRA' ...
    ' spectra_rob_eej p_st std_err st en periods wavenum n_data n_days']);
end;
%% Plot monthly data

days_in_prev_months = [0 31 59 90 120 151 181 212 243 273 304 334 366];

for i = 1:12,
st = days_in_prev_months(i);
en = days_in_prev_months(i+1);

    
    eval(['load /Users/manojnair/projects/tides/plots/' sprintf('DOY_%d_%d_LT_%d_%d_T_%s_%s',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29)) '_ROB_SPECTRA' ...
    ' spectra_rob_eej p_st std_err st en periods wavenum n_data n_days']);

spectra_rob_eej_mag = sqrt(spectra_rob_eej(:,:,2).^2+spectra_rob_eej(:,:,3).^2)*1e3;
spectra_rob_eej_phase = 24*((180/pi*angle(squeeze(complex(spectra_rob_eej(:,:,2),spectra_rob_eej(:,:,3)))))/360);
L = p_st(:,:,2) < p_test_lim & p_st(:,:,3) < p_test_lim; % Get the estimates for which both the real and imagianry parts are statistically significat
spectra_rob_eej_phase(~L) = 0; % Assign 0 to insigificant grid points
spectra_rob_eej_mag(~L) = 0;

spectra_monthly(i,:,:) = spectra_rob_eej_mag;
phase_monthly(i,:,:) = spectra_rob_eej_phase;

end;

% plot nonmigrating dirunal

fig=figure(1);
set(fig,'Position', [670         825        1007         265]);
[y2,ia2,ib2] = intersect(wavenum,[-4,-3,-2,-1,0,2,3,4]);
imagesc(datenum(2000,1:12,15),1:length(ia2), (squeeze(spectra_monthly(:,22,ia2)))');
set(gca,'YTick',1:length(ia2));
set(gca,'YTickLabel',['-4';'-3';'-2';'-1';' 0';' 2';' 3';' 4']);
title(['Non migrating dirunal' sprintf(' LT %d-%d p-limit %4.2f year %s-%s',lt_min,lt_max,p_test_lim,datestr(date_min,10),datestr(date_max,10))]);
ylabel('wavenumber');
h = colorbar; ylabel(h, 'EEJ amplitude mA/m');
datetick('x','mmm','keeplimits');
set(fig, 'PaperPositionMode', 'manual');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Dirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Amplitude'],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Dirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Amplitude'],'png');
close all;
% plot nonmigrating semi-dirunal
fig=figure(2);
set(fig,'Position', [670         825        1007         265]);

[y1,ia1,ib1] = intersect(periods,12);
[y2,ia2,ib2] = intersect(wavenum,[-3,-2,-1,0,1,3,4,5]);

imagesc(datenum(2000,1:12,15),1:length(ia2), (squeeze(spectra_monthly(:,ia1,ia2)))');
set(gca,'YTick',1:length(ia2));
%set(gca,'YTickLabel',['-4';'-3';'-2';'-1';' 0';' 1';' 3';' 4']);
set(gca,'YTickLabel',['-3';'-2';'-1';' 0';' 1';' 3';' 4';' 5']);
title(['Non migrating semi-dirunal ' sprintf(' LT %d-%d p-limit %4.2f year %s-%s',lt_min,lt_max,p_test_lim,datestr(date_min,10),datestr(date_max,10))]);
ylabel('wavenumber');
h = colorbar; ylabel(h, 'EEJ amplitude mA/m');
datetick('x','mmm','keeplimits');
set(fig, 'PaperPositionMode', 'manual');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Semidirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Amplitude'],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Semidirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Amplitude'],'png');

close all;
% plot nonmigrating teri-dirunal
fig=figure(3);
set(fig,'Position', [670         825        1007         265]);
[y1,ia1,ib1] = intersect(periods,8);

[y2,ia2,ib2] = intersect(wavenum,[0,1,2,4,5,6,7,8]);
imagesc(datenum(2000,1:12,15),1:length(ia2), (squeeze(spectra_monthly(:,ia1,ia2)))');
set(gca,'YTick',1:length(ia2));
%set(gca,'YTickLabel',['-4';'-3';'-2';'-1';' 0';' 1';' 2';' 4']);
set(gca,'YTickLabel',[' 0';' 1';' 2';' 3';' 4';' 5';' 6';' 7';' 8']);
title(['Non migrating teri-dirunal' sprintf(' LT %d-%d p-limit %4.2f year %s-%s',lt_min,lt_max,p_test_lim,datestr(date_min,10),datestr(date_max,10))]);
ylabel('wavenumber');
h = colorbar; ylabel(h, 'EEJ amplitude mA/m');
datetick('x','mmm','keeplimits');
set(fig, 'PaperPositionMode', 'manual');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Teridirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Amplitude'],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Teridirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Amplitude'],'png');
close all;

% plot nonmigrating q-dirunal
fig=figure(4);
set(fig,'Position', [670         825        1007         265]);

[y1,ia1,ib1] = intersect(periods,6);

[y2,ia2,ib2] = intersect(wavenum,[0,1,2,3,4,5,6,7,8]);
imagesc(datenum(2000,1:12,15),1:length(ia2), (squeeze(spectra_monthly(:,ia1,ia2)))');
set(gca,'YTick',1:length(ia2));
%set(gca,'YTickLabel',['-4';'-3';'-2';'-1';' 0';' 1';' 2';' 3']);
set(gca,'YTickLabel',[' 0';' 1';' 2';' 3';' 4';' 5';' 6';' 7';' 8']);
title(['Non migrating q-dirunal ' sprintf(' LT %d-%d p-limit %4.2f year %s-%s',lt_min,lt_max,p_test_lim,datestr(date_min,10),datestr(date_max,10))]);
ylabel('wavenumber');
h = colorbar; ylabel(h, 'EEJ amplitude mA/m');
datetick('x','mmm','keeplimits');
set(fig, 'PaperPositionMode', 'manual');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Qdirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Amplitude'],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Qdirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Amplitude'],'png');



% plot phase of the non-migrating tides

close all;
% plot nonmigrating dirunal
fig=figure(1);
set(fig,'Position', [670         825        1007         265]);
[y2,ia2,ib2] = intersect(wavenum,[-4,-3,-2,-1,0,2,3,4]);
imagesc(datenum(2000,1:12,15),1:length(ia2), squeeze(phase_monthly(:,22,ia2))');
set(gca,'YTick',1:length(ia2));
set(gca,'YTickLabel',['-4';'-3';'-2';'-1';' 0';' 2';' 3';' 4']);
title(['Non migrating dirunal ' sprintf(' LT %d-%d p-limit %4.2f year %s-%s',lt_min,lt_max,p_test_lim,datestr(date_min,10),datestr(date_max,10))]);
ylabel('wavenumber');
h = colorbar; ylabel(h, 'Phase hours');
datetick('x','mmm','keeplimits');
set(fig, 'PaperPositionMode', 'manual');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Dirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Phase'],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Dirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Phase'],'png');

close all;
% plot nonmigrating semi-dirunal
fig=figure(2);
set(fig,'Position', [670         825        1007         265]);

[y1,ia1,ib1] = intersect(periods,12);

[y2,ia2,ib2] = intersect(wavenum,[-3,-2,-1,0,1,3,4,5]);
imagesc(datenum(2000,1:12,15),1:length(ia2), (squeeze(phase_monthly(:,ia1,ia2)))');
set(gca,'YTick',1:length(ia2));
%set(gca,'YTickLabel',['-4';'-3';'-2';'-1';' 0';' 1';' 3';' 4']);
set(gca,'YTickLabel',['-3';'-2';'-1';' 0';' 1';' 3';' 4';' 5']);
title(['Non migrating semi-dirunal ' sprintf(' LT %d-%d p-limit %4.2f year %s-%s',lt_min,lt_max,p_test_lim,datestr(date_min,10),datestr(date_max,10))]);
ylabel('wavenumber');
h = colorbar; ylabel(h, 'Phase hours');
datetick('x','mmm','keeplimits');
set(fig, 'PaperPositionMode', 'manual');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Semidirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Phase'],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Semidirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Phase'],'png');

close all;
% plot nonmigrating teri-dirunal
fig=figure(3);
set(fig,'Position', [670         825        1007         265]);
[y1,ia1,ib1] = intersect(periods,8);

[y2,ia2,ib2] = intersect(wavenum,[0,1,2,4,5,6,7,8]);
imagesc(datenum(2000,1:12,15),1:length(ia2), (squeeze(phase_monthly(:,ia1,ia2)))');
set(gca,'YTick',1:length(ia2));
set(gca,'YTickLabel',[' 0';' 1';' 2';' 3';' 4';' 5';' 6';' 7';' 8']);
title(['Non migrating teri-dirunal ' sprintf(' LT %d-%d p-limit %4.2f year %s-%s',lt_min,lt_max,p_test_lim,datestr(date_min,10),datestr(date_max,10))]);
ylabel('wavenumber');
h = colorbar; ylabel(h, 'Phase hours');
datetick('x','mmm','keeplimits');
set(fig, 'PaperPositionMode', 'manual');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Teridirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Phase'],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Teridirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Phase'],'png');

close all;
% plot nonmigrating q-dirunal
fig=figure(4);
set(fig,'Position', [670         825        1007         265]);
[y1,ia1,ib1] = intersect(periods,6);
[y2,ia2,ib2] = intersect(wavenum,[0,1,2,4,5,6,7,8]);
imagesc(datenum(2000,1:12,15),1:length(ia2), (squeeze(phase_monthly(:,ia1,ia2)))');
set(gca,'YTick',1:length(ia2));
%set(gca,'YTickLabel',['-4';'-3';'-2';'-1';' 0';' 1';' 2';' 3']);
set(gca,'YTickLabel',[' 0';' 1';' 2';' 3';' 4';' 5';' 6';' 7';' 8']);
title(['Non migrating Q-dirunal ' sprintf(' LT %d-%d p-limit %4.2f year %s-%s',lt_min,lt_max,p_test_lim,datestr(date_min,10),datestr(date_max,10))]);
ylabel('wavenumber');
h = colorbar; ylabel(h, 'Phase hours');
datetick('x','mmm','keeplimits');
set(fig, 'PaperPositionMode', 'manual');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Qdirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Phase'],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Qdirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Phase'],'png');


close all;


%% Plot migrating tides using spectra_monthly and phase_monthly arrays from the above cell


[y1,ia1,ib1] = intersect(periods,[4,6,8,12,24]);
[y2,ia2,ib2] = intersect(wavenum,[1,2,3,4,5]);


j = 1;
for i = length(ia1):-1:1,
    
mig_mat(j,:) = (squeeze(spectra_monthly(:,ia1(i),ia2(j))));
mig_mat_phase(j,:) = (squeeze(phase_monthly(:,ia1(i),ia2(j))));
j = j+1;
end;

fig=figure(1);
set(fig,'Position', [670         825        1007         265]);
imagesc(datenum(2000,1:12,15),1:5, log10(mig_mat));
set(gca,'YTick',1:5);
set(gca,'YTickLabel',['1,1';'2,2';'3,3';'4,4';'5,5']);
title(['Amplitude of migrating waves ' sprintf(' LT %d-%d p-limit %4.2f year %s-%s',lt_min,lt_max,p_test_lim,datestr(date_min,10),datestr(date_max,10))] );
ylabel('n,s');
h = colorbar; ylabel(h, 'EEJ amplitude mA/m');
datetick('x','mmm','keeplimits');
set(fig, 'PaperPositionMode', 'manual');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Migrating_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Amplitude'],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Migrating_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Amplitude'],'png');


close all;
% phase of the tides
fig=figure(2);
set(fig,'Position', [670         825        1007         265]);
imagesc(datenum(2000,1:12,15),1:5, (mig_mat_phase));
set(gca,'YTick',1:5);
set(gca,'YTickLabel',['1,1';'2,2';'3,3';'4,4';'5,5']);
title(['Phase of migrating waves ' sprintf(' LT %d-%d p-limit %4.2f year %s-%s',lt_min,lt_max,p_test_lim,datestr(date_min,10),datestr(date_max,10))] );
ylabel('n,s');
h = colorbar; ylabel(h, 'Phase hours');
datetick('x','mmm','keeplimits');
set(fig, 'PaperPositionMode', 'manual');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Migrating_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Phase'],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Migrating_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Phase'],'png');

close all;


%% Plot the seasons data
season_st = [dayofyear(2000,05,01) dayofyear(2000,06,15) dayofyear(2000,11,01)];
season_en = [dayofyear(2000,08,30) dayofyear(2000,10,15) dayofyear(2000,02,29)];
fig = figure(1);
for seasons = 1:3,

 st = season_st(seasons);
 en = season_en(seasons);


 eval(['load /Users/manojnair/projects/tides/plots/' sprintf('DOY_%d_%d_LT_%d_%d_T_%s_%s',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29)) '_ROB_SPECTRA' ...
    ' spectra_rob_eej p_st std_err st en periods wavenum n_data n_days']);

spectra_rob_eej_mag = sqrt(spectra_rob_eej(:,:,2).^2+spectra_rob_eej(:,:,3).^2)*1e3;
spectra_rob_eej_phase = 180/pi*angle(squeeze(complex(spectra_rob_eej(:,:,2),spectra_rob_eej(:,:,3))));
L = p_st(:,:,2) < 0.05 & p_st(:,:,3) < p_test_lim; % Get the estimates for which both the real and imagianry parts are statistically significat
spectra_rob_eej_phase(~L) = 0; % Assign 0 to insigificant grid points
spectra_rob_eej_mag(~L) = 0;
subplot(length(season_st),2,(seasons-1)*2 + 1);
imagesc( wavenum, periods,log10(spectra_rob_eej_mag));
set(gca,'FontSize',10);
title( sprintf('DOY %d - %d,LT %d - %d,p-limit %4.2f,year %s-%s',st,en, lt_min,lt_max, p_test_lim,datestr(date_min,10),datestr(date_max,10)));
axis([-6 10 3 26]);
caxis([0,2]);
h = colorbar; ylabel(h, 'EEJ amplitude log(mA/m)');;
subplot(length(season_st),2,(seasons-1)*2 + 2);
imagesc( wavenum, periods,24 * ( spectra_rob_eej_phase./360 ));
axis([-6 10 3 26]);
caxis([-12,12]);
set(gca,'FontSize',10);
 h = colorbar; ylabel(h, 'Phase hours');
title( sprintf('DOY %d - %d,LT %d - %d,p-limit %4.2f,year %s-%s',st,en, lt_min,lt_max, p_test_lim,datestr(date_min,10),datestr(date_max,10)));

end;
set(fig, 'PaperPositionMode', 'manual');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Seasons_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10))],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Seasons_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10))],'png');

close all;
%% plot ALL data

st = 1;
en = 366;

 eval(['load /Users/manojnair/projects/tides/plots/' sprintf('DOY_%d_%d_LT_%d_%d_T_%s_%s',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29)) '_ROB_SPECTRA' ...
    ' spectra_rob_eej p_st std_err st en periods wavenum n_data n_days']);

spectra_rob_eej_mag = sqrt(spectra_rob_eej(:,:,2).^2+spectra_rob_eej(:,:,3).^2)*1e3;
spectra_rob_eej_phase = 180/pi*angle(squeeze(complex(spectra_rob_eej(:,:,2),spectra_rob_eej(:,:,3))));
L = p_st(:,:,2) < 0.05 & p_st(:,:,3) < p_test_lim; % Get the estimates for which both the real and imagianry parts are statistically significat
spectra_rob_eej_phase(~L) = 0; % Assign 0 to insigificant grid points
spectra_rob_eej_mag(~L) = 0;
figure(1);
imagesc( wavenum, periods,log10(spectra_rob_eej_mag));
set(gca,'FontSize',16);
title( sprintf('DOY %d - %d,LT %d - %d,p-limit %4.2f,year %s-%s',st,en, lt_min,lt_max, p_test_lim,datestr(date_min,10),datestr(date_max,10)));
axis([-6 10 3 26]);
caxis([0,2]);
h = colorbar; ylabel(h, 'EEJ amplitude log(mA/m)');
xlabel('wavenumber');
ylabel('period (hours)');

saveas(gcf, ['/Users/manojnair/projects/tides/plots/ALL_Amplitude_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10))],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/ALL_Amplitude_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10))],'png');

figure(2);
imagesc( wavenum, periods,24 * ( spectra_rob_eej_phase./360 ));
axis([-6 10 3 26]);
caxis([-12,12]);
set(gca,'FontSize',16);
h = colorbar; ylabel(h, 'Phase in hours');
title( sprintf('DOY %d - %d,LT %d - %d,p-limit %4.2f,year %s-%s',st,en, lt_min,lt_max, p_test_lim,datestr(date_min,10),datestr(date_max,10)));
set(fig, 'PaperPositionMode', 'manual');
xlabel('wavenumber');
ylabel('period (hours)');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/ALL_Phase_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10))],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/ALL_Phase_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10))],'png');

close all;

%% Print the results as a table
p_test_lim = 0.05; % Just to make sure that we don't print garbage !
fid = fopen('/Users/manojnair/projects/tides/plots/ionospheric_results.txt','at');

season_st = [dayofyear(2000,05,01) dayofyear(2000,06,15) dayofyear(2000,11,01)];
season_en = [dayofyear(2000,08,30) dayofyear(2000,10,15) dayofyear(2000,02,29)];

for seasons = 1:3,

 st = season_st(seasons);
 en = season_en(seasons);
 
fprintf(fid,'DOY_%d_%d_LT_%d_%d_T_%s_%s\n',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29));
eval(['load /Users/manojnair/projects/tides/plots/' sprintf('DOY_%d_%d_LT_%d_%d_T_%s_%s',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29)) '_ROB_SPECTRA' ...
    ' spectra_rob_eej p_st std_err st en periods wavenum n_data n_days']);
spectra_rob_eej_mag = sqrt(spectra_rob_eej(:,:,2).^2+spectra_rob_eej(:,:,3).^2)*1e3;
spectra_rob_eej_phase = 24*((180/pi*angle(squeeze(complex(spectra_rob_eej(:,:,2),spectra_rob_eej(:,:,3)))))/360);
L = p_st(:,:,2) < 0.05 & p_st(:,:,3) < p_test_lim; % Get the estimates for which both the real and imagianry parts are statistically significat
spectra_rob_eej_phase(~L) = 0; % Assign 0 to insigificant grid points
spectra_rob_eej_mag(~L) = 0;

[y1,ia1,ib1] = intersect(periods,[24,12,8,6,4]);
[y2,ia2,ib2] = intersect(wavenum,[-7:7]);

for i = 1:length(ia1),
    for j = 1:length(ia2),
        
        if p_st(ia1(i),ia2(j),2) <0.05 & p_st(ia1(i),ia2(j),3) <0.05,
            fprintf(fid,'%2d,%2d   %6.2f %6.2f %6.2f %6.2f %6.2f %5.2f\n', ...
                periods(ia1(i)),wavenum(ia2(j)), ...
                spectra_rob_eej(ia1(i),ia2(j),2)*1e3, ...
                spectra_rob_eej(ia1(i),ia2(j),3)*1e3, ...
                std_err(ia1(i),ia2(j),2)*1e3, ...
                std_err(ia1(i),ia2(j),3)*1e3, ...
                spectra_rob_eej_mag(ia1(i),ia2(j)), ...
                spectra_rob_eej_phase(ia1(i),ia2(j)));
        end;
    end;
end;
                
    
end;


fclose(fid);


%% print the results for all data
p_test_lim = 0.05; % Just to make sure that we don't print garbage !
fid = fopen('/Users/manojnair/projects/tides/plots/ionospheric_results.txt','at');
st = 1;
en = 366;

 eval(['load /Users/manojnair/projects/tides/plots/' sprintf('DOY_%d_%d_LT_%d_%d_T_%s_%s',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29)) '_ROB_SPECTRA' ...
    ' spectra_rob_eej p_st std_err st en periods wavenum n_data n_days']);
fprintf(fid,'DOY_%d_%d_LT_%d_%d_T_%s_%s\n',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29));
spectra_rob_eej_mag = sqrt(spectra_rob_eej(:,:,2).^2+spectra_rob_eej(:,:,3).^2)*1e3;
spectra_rob_eej_phase = 24*((180/pi*angle(squeeze(complex(spectra_rob_eej(:,:,2),spectra_rob_eej(:,:,3)))))/360);
L = p_st(:,:,2) < 0.05 & p_st(:,:,3) < 0.05; % Get the estimates for which both the real and imagianry parts are statistically significat
spectra_rob_eej_phase(~L) = 0; % Assign 0 to insigificant grid points
spectra_rob_eej_mag(~L) = 0;

[y1,ia1,ib1] = intersect(periods,[24,12,8,6,4]);
[y2,ia2,ib2] = intersect(wavenum,[-7:7]);

for i = 1:length(ia1),
    for j = 1:length(ia2),
        
        if p_st(ia1(i),ia2(j),2) <0.05 & p_st(ia1(i),ia2(j),3) <0.05,
            fprintf(fid,'%2d,%2d   %6.2f %6.2f %6.2f %6.2f %6.2f %5.2f\n', ...
                periods(ia1(i)),wavenum(ia2(j)), ...
                spectra_rob_eej(ia1(i),ia2(j),2)*1e3, ...
                spectra_rob_eej(ia1(i),ia2(j),3)*1e3, ...
                std_err(ia1(i),ia2(j),2)*1e3, ...
                std_err(ia1(i),ia2(j),3)*1e3, ...
                spectra_rob_eej_mag(ia1(i),ia2(j)), ...
                spectra_rob_eej_phase(ia1(i),ia2(j)));
        end;
    end;
end;
fclose(fid);


% 
% %%
% %testing whether ordered longitude is necessary
% 
% [y,ia] = sort(eef(:,2));
% 
% x_lon = [ cos(2*pi*((lonnorm(ia)'*2))) ...
%           sin(2*pi*((lonnorm(ia)'*2)))];  
% 
% spectra_ln = x_lon\eef(ia,6);
% 
% % now invert without ordering
% 
% x_lon_n = [ cos(2*pi*((lonnorm'*3.5))) ...
%           sin(2*pi*((lonnorm'*3.5)))];  
% 
% spectra_ln_n = x_lon_n\eef(:,6);
% 
% %inference. There is no difference between ordered or non-ordered inversion
% 
% % finding out the wave structure of M2  tides
% 
% 
% periods = 12.415:0.001:12.43 ; 
% wavenum = -8:1:8;
% clear spectra_ls spectra_rob_eej spectra_ls_eej
% for i = 1:length(periods),
%     for j = 1 : length(wavenum),
%     
% prd = periods(i);
% 
% x_m2 = [ cos(2*pi*((eef(:,1)/( prd/24)) + (lonnorm'*wavenum(j)))) ...
%          sin(2*pi*((eef(:,1)/( prd/24)) + (lonnorm'*wavenum(j))))];  
% [spectra_rob_eej(i,j,:), stats(i,j) ] = robustfit(x_m2,eef(:,10));
% %spectra_ls_eej(i,j,:) = x_m2\(eef(LL,10));
% end;
% end;
% 

%% Joint CHAMP EEF and Jicamarca ISR EEF analysis
% The idea is to include the Icamarca ISR data to the inversion. This
% should remove the "nigh-time" ambiguity in the wave fitting. Inference.
% By including the ISR data, the wavemodes are no longer discrete, but
% continous signals, with the highest at the migrating tides mode. For
% example, the relative differences between DE3 and SE2 is not very clear
% now. However, one should also look whether Jicamarca data is overwhelming the
% inversion process. For example, when I limited jicamarca data to just
% 5000 samples, the spectral image became much similar to the one I got
% without including Jicamarca data.

load /data/backup/mnair/ace/jcamarca_isr_fejer.mat eef;
jeef = eef;
load /data/backup/mnair/longp/eef_data_2000_2009.mat eef;


jeef(:,1) = jeef(:,1) + datenum(2000,1,1);

% running mean of Jicamarca ISR 

jic_eef = runmean(jeef(1:5000,6),1); % 


% add jicamarca data to champ data

t = [eef(:,1)' jeef(1:5000,1)']';
y = [eef(:,6)' jic_eef']';
l = [eef(:,2)' ones([1,length(jic_eef)])*-76.8];
lt = [eef(:,4)' jeef(1:5000,4)'];

lt_min = 7;
lt_max = 17;
date_min = datenum(2000,08,1);
date_max = datenum(2005,07,31);
p_test_lim = 0.05;
periods = 3:1:30;
wavenum = -10:1:10;


LL =  t > date_min & t < date_max;

[n_data,n_days]=hist(doyvec(LL));

% make normalized latitudes between 0 (= 0) and 360 (=1)

if exist('lonnorm'),
    clear lonnorm;
end;

selected_lon = l(LL);
L = selected_lon < 0;
lonnorm(L) = 360 + selected_lon(L) ;
lonnorm(~L) =  selected_lon(~L) ;
lonnorm = lonnorm/360;


clear spectra_ls spectra_rob_eej spectra_ls_eej p_st std_err stats
for i = 1:length(periods),
    for j = 1 : length(wavenum),
    
prd = periods(i);

x_m2 = [ cos(2*pi*((t(LL)/( prd/24)) + (lonnorm'*wavenum(j)))) ...
         sin(2*pi*((t(LL)/( prd/24)) + (lonnorm'*wavenum(j))))];  
%[spectra_rob_eej_jic(i,j,:), stats(i,j) ] = robustfit(x_m2,y(LL));
% p_st(i,j,:) = stats(i,j).p;
% std_err(i,j,:) = stats(i,j).se;
spectra_ls_eej_jic2(i,j,:) = x_m2\(y(LL));
end;
end;
